#include "pdf.h"
#include "tensor.h"
#include "vector_tools.h"
#include <stdexcept>
#include <fstream>
#include <sstream>
#include <vector>

//return sign of double
inline int sign(double input)
{
  if(input>0) return 1;
  else return -1;
}

//core function to do the interpolation
double PDF::operator()(const vector<double>& ntuple, double* err) const
{
  //make sure input has the same dimension
  assert(ntuple.size()==data.rank);

  //grid coordinate = relative coordinate with respect
  //to the bin in which the extrapolation is carried out
  vector<double> grid_coordinate(data.rank);
  vector<int> grid_index(data.rank);

  //DEBUG
  /*
  cout<<"PDF: input ntuple "<<ntuple<<endl;
  */

  //see which direction to compute derivative
  //+ or -
  vector<int> bin_direction(data.rank);

  //find the grid coordinate
  for(unsigned int i=0; i<ntuple.size(); ++i)
    {
      //compute bin_coordinate
      double temp_coordinate=(ntuple[i]-min[i])*one_over_width[i];

      //compute which bin this falls in
      int temp_bin=static_cast<int>(temp_coordinate);

      //shift grid coordinate to -0.5
      //since the first bin stores the value at bin coordinate 0.5
      //not 0
      temp_coordinate-=0.5;
      
      //minimum is -0.5, 0 is the middle of the first bin
      if(temp_coordinate < -0.5) temp_coordinate=-0.5;

      //maximum is nbins[i]-0.5, since nbins[i]-1 
      //stores values not at the maximum but middle of last bin
      else if(temp_coordinate > nbins[i]-0.5) 
	temp_coordinate=nbins[i]-0.5;

      //minimum bin is 0
      if(temp_bin<0) temp_bin=0;

      //maximum bin is nbins[i]-1
      else if(temp_bin>=nbins[i])
	temp_bin=nbins[i]-1;

      //get the bin direction
      //this is useful mainly for points near the boundary
      int temp_bin_dir=sign(temp_coordinate-temp_bin);

      //if bin is at lower boundary
      //and evaluation is done below the middle of the boundary bound
      //this will demand special treatment
      //set direction to zero for now
      //since there is no more data below
      if(temp_bin==0 && temp_bin_dir<0)
	temp_bin_dir=0;	

      //same thing for maximum
      //there is no more bin above
      else if(temp_bin==nbins[i]-1 && temp_bin_dir>0)
	temp_bin_dir=0;

      //now populate the vectors
      grid_coordinate[i]=temp_coordinate-temp_bin;
      grid_index[i]=temp_bin;
      bin_direction[i]=temp_bin_dir;
    }

  //DEBUG
  /*
  cout<<"PDF: grid_coordinate "<<grid_coordinate<<endl;
  cout<<"PDF: grid_index "<<grid_index<<endl;
  cout<<"PDF: bin_direction "<<bin_direction<<endl;
  */
  
  //now we know which bin to look for info and
  //what bin to look for derivative
  //fist compute the gradient
  vector<double> gradient(data.rank);

  //do the same for err
  vector<double> gradient_err;

  //only initialize err if it is present
  if(err)
    gradient_err.resize(data.rank);

  double bin_value=data[grid_index];
  double err_value=data_err[grid_index];

  for(int i=0; i<gradient.size(); ++i)
    {
      vector<int> current_index=grid_index;

      //special treatment for boundary
      if(bin_direction[i]==0)
	{
	  bin_direction[i]=sign(grid_coordinate[i]);
	  //simply extrapolate gradient assuming zero at boundary
	  gradient[i]=bin_direction[i]*(0.0-bin_value);

	  if(err)
	    gradient_err[i]=bin_direction[i]*(0.0-err_value);
	}

      //for other cases
      //gradient is easily calculated
      else
	{
	//in the ith direction
	//we want to look in either the next
	//or previous cell
	current_index[i]+=bin_direction[i];
	
	//Tensor class takes input indices
	//and grab their values      
	gradient[i]= bin_direction[i]*(data[current_index]-bin_value);

	if(err)
	  gradient_err[i]= bin_direction[i]*(data_err[current_index]-err_value);
	}
    }

  //DEBUG
  /*
  cout<<"PDF: gradient "<<gradient<<endl;
  cout<<"PDF: bin value "<<bin_value<<endl;
  */

  //store error if pointer is specified
  //no need to do complicated interpolation for error
  double result=(gradient*grid_coordinate) + bin_value;

  if(err)
    {
      //do a little linear interpolation
      (*err)=(gradient_err*grid_coordinate) + err_value;

      //now if err is very large (low statistics)
      //or if pdf value is extremely low
      //we demand it has to be at least 
      //larger than min error
      /*
      if( (*err)<min_err && 
	  (result < 1e-21 || data_err[grid_index] > 0.5 ))
	(*err)=min_err;
      */
    }
  //now gradient is computed
  //we can interpolate!
  //slick use of vector! dotproduct
  return result;
}

//fully evaluate PDF and gradient
double PDF::FullEvaluate(const vector<double>& ntuple, double* err,
			 vector<double>& grad_output) const
{
  //make sure input has the same dimension
  assert(ntuple.size()==data.rank);
  assert(grad_output.size()==data.rank);

  //grid coordinate = relative coordinate with respect
  //to the bin in which the extrapolation is carried out
  vector<double> grid_coordinate(data.rank);
  vector<int> grid_index(data.rank);

  //DEBUG
  /*
  cout<<"PDF: input ntuple "<<ntuple<<endl;
  */

  //see which direction to compute derivative
  //+ or -
  vector<int> bin_direction(data.rank);

  //find the grid coordinate
  for(unsigned int i=0; i<ntuple.size(); ++i)
    {
      //compute bin_coordinate
      double temp_coordinate=(ntuple[i]-min[i])*one_over_width[i];

      //compute which bin this falls in
      int temp_bin=static_cast<int>(temp_coordinate);

      //shift grid coordinate to -0.5
      //since the first bin stores the value at bin coordinate 0.5
      //not 0
      temp_coordinate-=0.5;
      
      //minimum is -0.5, 0 is the middle of the first bin
      if(temp_coordinate < -0.5) temp_coordinate=-0.5;

      //maximum is nbins[i]-0.5, since nbins[i]-1 
      //stores values not at the maximum but middle of last bin
      else if(temp_coordinate > nbins[i]-0.5) 
	temp_coordinate=nbins[i]-0.5;

      //minimum bin is 0
      if(temp_bin<0) temp_bin=0;

      //maximum bin is nbins[i]-1
      else if(temp_bin>=nbins[i])
	temp_bin=nbins[i]-1;

      //get the bin direction
      //this is useful mainly for points near the boundary
      int temp_bin_dir=sign(temp_coordinate-temp_bin);

      //if bin is at lower boundary
      //and evaluation is done below the middle of the boundary bound
      //this will demand special treatment
      //set direction to zero for now
      //since there is no more data below
      if(temp_bin==0 && temp_bin_dir<0)
	temp_bin_dir=0;	

      //same thing for maximum
      //there is no more bin above
      else if(temp_bin==nbins[i]-1 && temp_bin_dir>0)
	temp_bin_dir=0;

      //now populate the vectors
      grid_coordinate[i]=temp_coordinate-temp_bin;
      grid_index[i]=temp_bin;
      bin_direction[i]=temp_bin_dir;
    }

  //DEBUG
  /*
  cout<<"PDF: grid_coordinate "<<grid_coordinate<<endl;
  cout<<"PDF: grid_index "<<grid_index<<endl;
  cout<<"PDF: bin_direction "<<bin_direction<<endl;
  */
  
  //now we know which bin to look for info and
  //what bin to look for derivative
  //fist compute the gradient
  vector<double> gradient(data.rank);
  double bin_value=data[grid_index];

  for(int i=0; i<gradient.size(); ++i)
    {
      vector<int> current_index=grid_index;

      //special treatment for boundary
      //extrapolate from zero to value
      /*
      if(bin_direction[i]==0)
	{
	  bin_direction[i]=sign(grid_coordinate[i]);
	  //simply extrapolate gradient assuming zero at boundary
	  gradient[i]=bin_direction[i]*(0.0-bin_value);
	}
      */


      //for other cases
      //gradient is easily calculated
	//in the ith direction
	//we want to look in either the next
	//or previous cell
	current_index[i]+=bin_direction[i];
	
	//Tensor class takes input indices
	//and grab their values      
	gradient[i]= bin_direction[i]*(data[current_index]-bin_value);
    }

  //DEBUG
  /*
  cout<<"PDF: gradient "<<gradient<<endl;
  cout<<"PDF: bin value "<<bin_value<<endl;
  */

  //store the gradient in grad_output
  grad_output=prod(gradient, one_over_width);

  //store error if pointer is specified
  //no need to do complicated interpolation for error
  double result=(gradient*grid_coordinate) + bin_value;

  if(err!=NULL)
    {
      (*err)=data_err[grid_index];

      //now if err is very large (low statistics)
      //or if pdf value is extremely low
      //we demand it has to be at least 
      //larger than min error
      if( (*err)<min_err && 
	  (result < 1e-21 || data_err[grid_index] > 0.5 ))
	(*err)=min_err;
    }
  //now gradient is computed
  //we can interpolate!
  //slick use of vector! dotproduct
  return result;
}

//function to output PDF to file
bool PDF::Write(const char* output) const
{
  //open a file
  ofstream fout(output);

  if(!fout.good())
    return false;
  

  //first output some formatting
  fout<<"<num>"<<endl;
  
  //then cout numbers per line
  //by default it is 10
  int n_per_line=10;

  fout<<n_per_line<<endl;
  fout<<"</num>"<<endl<<endl;

  //output rank
  fout<<"<rank>"<<endl;
  fout<<Rank()<<endl<<"</rank>"<<endl<<endl;

  //now cout number of bins
  fout<<"<nbins>"<<endl;

  for(unsigned int i=0; i<nbins.size(); ++i)
    fout<<nbins[i]<<" ";
  fout<<endl<<"</nbins>"<<endl<<endl;

  //now cout max and min
  fout<<"<max>"<<endl;

  for(unsigned int i=0; i<max.size(); ++i)
    fout<<max[i]<<" ";
  fout<<endl<<"</max>"<<endl<<endl;

  fout<<"<min>"<<endl;

  for(unsigned int i=0; i<min.size(); ++i)
    fout<<min[i]<<" ";
  fout<<endl<<"</min>"<<endl<<endl;
  
  //the pdf itself
  fout<<"<pdf>";
  for(unsigned int i=0; i<data.size(); ++i)
    {
      if(i% n_per_line ==0)
	fout<<endl;
      fout<<data[i]<<" ";
    }
  fout<<endl<<"</pdf>"<<endl<<endl;

  //bias
  fout<<"<bias>";
  for(unsigned int i=0; i<data_bias.size(); ++i)
    {
      if(i% n_per_line ==0)
	fout<<endl;
      fout<<data_bias[i]<<" ";
    }
  fout<<endl<<"</bias>"<<endl<<endl;

  //variance
  fout<<"<var>";
  for(unsigned int i=0; i<data_var.size(); ++i)
    {
      if(i% n_per_line ==0)
	fout<<endl;
      fout<<data_var[i]<<" ";
    }
  fout<<endl<<"</var>"<<endl<<endl;

  //err
  fout<<"<err>";
  for(unsigned int i=0; i<data_err.size(); ++i)
    {
      if(i% n_per_line ==0)
	fout<<endl;
      fout<<data_err[i]<<" ";
    }
  fout<<endl<<"</err>"<<endl<<endl;
  

  //min_err
  fout<<"<minerr>"<<endl
      <<min_err<<endl
      <<"</minerr>"<<endl<<endl;
  

  fout.close();
  return true;
  
}

//function to read PDF from file
bool PDF::Read(const char* input)
{
  //open file input
  ifstream fin(input);

  if(!fin.good())
    return false;
  
  string line;

  //see which block we are reading
  // 0 = nothing
  // 1 = num
  // 2 = nbins
  // 3 = max
  // 4 = min
  // 5 = pdf
  // 6 = bias
  // 7 = variance
  // 8 = err
  // 9 = minerr
  // 10 = rank
  
  int states=0;
  int n_per_line=0;
  int rank=0;

  //various indices
  int n_nbins=0;
  int n_max=0;
  int n_min=0;
  int n_pdf=0;
  int n_bias=0;
  int n_var=0;
  int n_err=0;
  
  //make sure all the variables are present
  bool bstates[10];

  while(fin.good())
    {
      //get the line
      getline(fin, line);

      if(line.length() < 1) continue;
      
      //if we are not inbetween blocks
      if(states == 0)
	{
	  //try to look for open bracket
	  if(line[0]=='<')
	    {
	      if(line=="<num>")
		{states=1;}
	      else if(line=="<nbins>")
		{
		  states=2;
		  nbins.resize(rank);
		}
	      else if(line=="<max>")
		{
		  states=3;
		  max.resize(rank);
		}
	      else if(line=="<min>")
		{
		  states=4;
		  min.resize(rank);
		}
	      else if(line=="<pdf>")
		{
		  states=5;
		  data=Tensor<double>(nbins);
		}
	      else if(line=="<bias>")
		{
		  states=6;
		  data_bias=Tensor<double>(nbins);
		}
	      else if(line=="<var>")
		{
		  states=7;
		  data_var=Tensor<double>(nbins);
		}
	      else if(line=="<err>")
		{
		  states=8;
		  data_err=Tensor<double>(nbins);
		}
	      else if(line=="<minerr>")
		{states=9;}
	      else if(line=="<rank>")
		{states=10;}
	      
	      //see if the entries per line is defined    
	      if(states > 1 && n_per_line==0)
		{
		  cerr<<"ERROR: number of entries per line ill-defined in PDF"
		      <<endl;
		  return false;
		}
	      if(states == 2 && rank==0)
		{
		  cerr<<"ERROR: cannot determined rank"
		      <<endl;
		  return false;
		}
	      
	    }
	  continue;
	}

      //now states must be > 0
      //if the block is ending
      if(line[0]=='<' && line.length() >2 && line[1]=='/')
	{
	  bstates[states-1]=true;
	  states=0;
	  continue;
	}
      
      stringstream ss(line);
      
      //if the block is not ending, grab numbers
      switch(states)
	{
	  //get entries per line
	case 1:
	  {
	    ss >> n_per_line;
	    if(n_per_line <= 0 || n_per_line > 100)
	      {
		cerr<<"ERROR: number of entries per line ill-defined in PDF"
		    <<endl;
		return false;
	      }
	    break;
	  }
	  
	  //get the rank of PDF
	case 10:
	  {
	    ss >> rank;
	    if(rank <= 0 || rank > 100) 
	      {
		cerr<<"ERROR: rank ill-defined in PDF "<<rank
		    <<endl;
		return false;
	      }
	    break;
	  }   
	  
	  //populate nbins
	case 2:
	  {
	    for(unsigned int i=0; i<n_per_line && n_nbins<rank; ++i)
	      {
		ss >> nbins[n_nbins];
		++n_nbins;
	      }
	    break;
	  }
	  
	  //populate max
	case 3:
	  {
	    for(unsigned int i=0; i<n_per_line && n_max<rank; ++i)
	      {
		ss >> max[n_max];
		++n_max;
	      }
	    break;
	  }
	  
	  //populate min
	case 4:
	  {
	    for(unsigned int i=0; i<n_per_line && n_min<rank; ++i)
	      {
		ss >> min[n_min];
		++n_min;
	      }
	    break;
	  }
	  
	  //populate pdf
	case 5:
	  {
	    for(unsigned int i=0; i<n_per_line && n_pdf<data.size(); ++i)
	      {
		ss >> data[n_pdf];
		++n_pdf;
	      }
	    break;
	  }
	  
	  //populate bias
	case 6:
	  {
	    for(unsigned int i=0; i<n_per_line && n_bias<data_bias.size();
		++i)
	      {
		ss >> data_bias[n_bias];
		++n_bias;
	      }
	    break;
	  }
	  
	  //populate variance
	case 7:
	  {
	    for(unsigned int i=0; i<n_per_line && n_var<data_var.size();
		++i)
	      {
		ss >> data_var[n_var];
		++n_var;
	      }
	    break;
	  }
	  
	  //populate error
	case 8:
	  {
	    for(unsigned int i=0; i<n_per_line && n_err<data_err.size();
		++i)
	      {
		ss >> data_err[n_err];
		++n_err;
	      }
	    break;
	  }
	  
	  //min error
	case 9:
	  {
	    ss >> min_err;
	    break;
	  }

	default:
	  {
	    //some crazy error?
	    cerr<<"ERROR: non-defined states in read PDF??"<<endl;
	    return false;
	  } 
	}//done switch
      
    }//done while loop

  
  //get the width
  width.resize(rank);
  one_over_width.resize(rank);
  for(unsigned int i=0; i<rank; ++i)
    {
      width[i] = (max[i] - min[i])/nbins[i];
      one_over_width[i] = (1.0/width[i]);
    }
  
  //see if any states are empty
  for(int i=0; i<10; i++)
    if(!bstates[i])
      {
	cerr<<"ERROR: not enough info in input PDF text file, states: "
	    <<i+1<<endl;
	return false;
      }
  
  //PDF is now valid!
  valid=true;
  fin.close();
  return true;
}


PDF PDF::Project(const vector<Projection> input) const
{
  //first we have to determine the rank and size of new PDF
  PDF newpdf;
  if(input.size()!=data.rank)
    throw std::runtime_error
      ("The input to PDF::Project is not the right size.");
  
  //first check the new rank
  vector<double> newmin;
  vector<double> newmax;
  vector<double> newprecision;
  
  for(unsigned int i=0; i<input.size(); ++i)
    {
      if(input[i]==PDF::KEEP)
	{
	  newmin.push_back(min[i]);
	  newmax.push_back(max[i]);	  
	  newprecision.push_back((max[i]-min[i])*one_over_width[i]);
	}
    }
  
  //now we have all the new dimensions
  //initialize new pdf
  newpdf.Initialize(newmin.size(), &newmin[0], &newmax[0], &newprecision[0]);
  
  //we need to know which indices in the tensor to sum over
  
  //keep track of how many cells need to be summed over
  int n_integrate=1;
  
  //index vector to convert from newpdf to oldpdf indices
  vector<int> new_to_old(newmin.size());
  
  //index vector to convert the indices for
  //the integrated parameters to indices in the oldpdf
  vector<int> integrate_to_old(data.rank-newmin.size());
  
  //also we'll need a tensor to represent the orthogonal space
  vector<int> orth_length(data.rank-newmin.size());
  
  //now figure out what these converters are
  int new_to_old_index=0;
  int integrate_to_old_index=0;
  
  for(unsigned int i=0; i<input.size(); ++i)
    {
      //if the variable is to be kept
      //get the convertion factor to new_to_old
      if(input[i]==PDF::KEEP)
	{
	  new_to_old[new_to_old_index]= data.index_helper[i];
	  ++new_to_old_index;
	}
      //if the variable is not kept
      //give it to integrate_to_old
      else
	{
	  integrate_to_old[integrate_to_old_index]= data.index_helper[i];
	  orth_length[integrate_to_old_index]=data.length[i];
	  ++integrate_to_old_index;
	  //we have more dimensions to integrate
	  n_integrate*=data.length[i];
	}
    }

  //a slight detour, after getting n_integrate
  //we can set the new minimum error for the projected PDF
  newpdf.min_err=min_err/(n_integrate);

  //we need a tensor representing the orthogonal space
  //just for the function IndexVec
  Tensor<double> orth(orth_length);

  //here we have to perform the integration
  //first loop over all the cells in data
  for(unsigned int i=0; i<newpdf.data.size(); ++i)
    {
      //1. first do a simple sum
      //2. then normalize to one again
      //3. and compute new errors

      //get an initial position in the oldpdf first
      int initial_pos = (new_to_old * newpdf.data.IndexVec(i));
      
      //now sum over the bins
      for(int j=0; j<n_integrate; ++j)
	{
	  //get the location in the oldpdf to add
	  int pos= initial_pos + (integrate_to_old * orth.IndexVec(j));
	  
	  //now increment the PDF
	  newpdf.data[i]+= data[pos];
	  
	  //increment bias (in percentage)
	  newpdf.data_bias[i]+= data_bias[pos]*data[pos];

	  //variance add in quadrature
	  newpdf.data_var[i]+= data_var[pos]*data[pos]*data[pos];
	}

      //normalize the error (get percentage)
      //zero we do nothing
      if(newpdf.data[i] == 0)
	{
	  newpdf.data_bias[i]=0;
	  newpdf.data_var[i]=0;
	  newpdf.data_err[i]=0;
	}
      //now compute total error
      else
	{
	  newpdf.data_bias[i]/= newpdf.data[i];
	  newpdf.data_var[i]/= (newpdf.data[i]*newpdf.data[i]);

	  newpdf.data_err[i]= 
	    sqrt(newpdf.data_var[i] + newpdf.data_bias[i]*newpdf.data_bias[i]);
	} 
    }

  //normalize pdf
  //get the volume of cell
  double volume=1.0;
  for(unsigned int i=0; i<newpdf.width.size(); ++i)
    volume*=newpdf.width[i];

  double norm= summation(newpdf.data.ary)*volume;
  
  //if norm is identically zero, we have a problem
  if(norm ==0)
    throw std::runtime_error
      ("When Projecting PDF, integral of pdf is identically zero.");

  //rescale the whole pdf
  newpdf.data.ary /= norm;

  //we now have a new PDF!
  newpdf.valid=true;
  return newpdf;
}

void PDF::Initialize(int nvar, const double* min_, const double* max_, 
		     const double* precision_, double min_err_)
{
  //make sure inputs are well defined
  assert(nvar>0);
  assert(min_!=NULL);
  assert(max_!=NULL);
  
  min.resize(nvar);
  max.resize(nvar);
  one_over_width.resize(nvar);
  width.resize(nvar);
  nbins.resize(nvar);
  
  //set minimum PDF error for low statistics regions
  min_err=min_err_;

  //size for covariance matrix and data
  int datasize=1;

  for(int i=0; i<nvar; i++)
    {
      //set min and max values
      min[i]=min_[i];
      max[i]=max_[i];

      //set the number of bins
      if(precision_==NULL)
	{
	  nbins[i]=ceil((max[i]-min[i])/0.01);
	}
      else
	{
	  assert(precision_[i]>0);
	  nbins[i]=ceil((max[i]-min[i])/precision_[i]);
	}
      
      //the size of each dimension must be greater than zero
      assert(nbins[i]>0);

      datasize*=nbins[i];
      //now we may compute the width
      one_over_width[i]=static_cast<double>(nbins[i])/(max[i]-min[i]);
      width[i]=(max[i]-min[i])/static_cast<double>(nbins[i]);
    }

  //initialize the data grid
  data=Tensor<double>(nbins);
  data_bias=Tensor<double>(nbins);
  data_var=Tensor<double>(nbins);
  data_err=Tensor<double>(nbins);
}


ostream& operator<<(ostream& os, const PDF& pdf)
{
  os<<endl<<"PDF info: (size="<<pdf.data.size()<<")"<<endl;
  for(int i=0; i<pdf.data.size(); ++i)
    os<<pdf.data[i]<<endl;

  return os;
}

//code for MC_Generator sub class in PDF
//initialize random number generator at run time
MTRand_int32 PDF::MC_Generator::irand(time(0));

PDF::MC_Generator::MC_Generator(const vector<MC_Flag>& flag_, const PDF* pdf_):
flag(flag_), pdf(pdf_)
{
  if(pdf==NULL)
    throw std::runtime_error
      ("PDF is null when creating MC_Generator");

  if(flag.size() != pdf->Rank())
    throw std::runtime_error
      ("MC_Generator input does not match rank of PDF");
  
  input_size=0;

  for(unsigned i=0; i<flag.size(); ++i)
    {
      if(flag[i]==PDF::INPUT)
	++input_size;
      else
	index.push_back(i);
    }

  if(input_size == pdf->Rank())
    throw std::runtime_error
      ("All MC_Generator flags are input, there is nothing to generate!");
}

//generate event
vector<PDF::MC_Event> PDF::MC_Generator::Generate(const vector<double>& input, 
					  int nevt) const
{

  //the number of input must match the input_size
  assert(input.size()==input_size);

  //make sure the PDF MC_Generator refers to exists
  assert(pdf!=NULL);

  //initialize result
  vector<PDF::MC_Event> result;
  result.reserve(nevt);

  //initialize an input and output vector
  vector<double> base_input(pdf->Rank());
  vector<double> base_output(pdf->Rank() - input_size);

  //populate input values

  //index for the input vector
  int input_index=0;

  //index to see if current i is an input index
  vector<int>::const_iterator current_index_check=index.begin();
  for(unsigned int i=0; i<base_input.size(); ++i)
    {
      //if the index i is an input, ignore
      if((*current_index_check)==i)
	{
	  if(current_index_check!=index.end())
	    ++current_index_check;
	  continue;
	}

      //else it's an input, simply grab it
      else
	{
	  base_input[i]=input[input_index];
	  //on to the next input
	  //potential problem if input_index goes out of bound??
	  ++input_index;
	}
    }

  //first loop over however many events
  for(int i=0; i<nevt; i++)
    {

      //loop over all the output that needs Monte Carlo
      for(unsigned j=0; j<index.size(); ++j)
	{
	  base_input[index[j]]=
	    pdf->min[index[j]] + 
	    (pdf->max[index[j]]-pdf->min[index[j]])*irand.rand_double();
	  
	  base_output[j]=base_input[index[j]];
	}      
      //evaluate the pdf at that point
      double output_err;
      double output_value=(*pdf)(base_input, &output_err);

      //we have weight, error, time to push_back the new event!
      result.push_back(PDF::MC_Event(base_output, output_err, output_value));
    }

  //now we have to rescale all the weight so they all sum to 1
  //first get the sum of all weight
  double sum_of_weight=0;
  for(unsigned int i=0; i<result.size(); ++i)
    sum_of_weight+=result[i].weight;

  //see if the sum of weights are zero (something wrong??)
  if(sum_of_weight==0)
    cout<<"ERROR: pdf_values from MC_Generator are all zero, results may not make sense"<<endl;
  
  if(sum_of_weight>0)
    {
      //get 1/sum first to make calculation more efficient
      sum_of_weight=1.0/sum_of_weight;
      
      for(unsigned int i=0; i<result.size(); ++i)
	{
	  result[i].weight*=sum_of_weight;
	  result[i].err*=sum_of_weight;
	}
    }

  //and that's it! return the list of events
  return result;
}



vector<PDF::MC_Event> PDF::MC_Generator::GenerateGrad
(const vector<double>& input, int nevt) const
{
  
  //the number of input must match the input_size
  assert(input.size()==input_size);

  //make sure the PDF MC_Generator refers to exists
  assert(pdf!=NULL);

  //initialize result
  vector<PDF::MC_Event> result;
  result.reserve(nevt);

  //initialize an input and output vector
  vector<double> base_input(pdf->Rank());
  vector<double> base_output(pdf->Rank() - input_size);

  //populate input values

  //index for the input vector
  int input_index=0;

  //gradient pointer
  vector<double> gradient(pdf->Rank());

  //we also need the volume for the output
  //used to normalize the gradient
  double gradient_volume=1.0;

  for(unsigned int i=0; i<base_input.size(); ++i)
    {
      if(flag[i]==INPUT)
	{
	  base_input[i]=input[input_index];
	  //on to the next input
	  //potential problem if input_index goes out of bound??
	  ++input_index;
	}
      else
	{
	  //get the gradient_volume
	  gradient_volume*=(pdf->max[i]-pdf->min[i]);
	}
    }

  //first loop over however many events
  for(int i=0; i<nevt; i++)
    {

      //loop over all the output that needs Monte Carlo
      for(unsigned j=0; j<index.size(); ++j)
	{
	  base_input[index[j]]=
	    pdf->min[index[j]] + 
	    (pdf->max[index[j]]-pdf->min[index[j]])*irand.rand_double();
	  
	  base_output[j]=base_input[index[j]];
	}      
      //evaluate the pdf at that point
      double output_err;
      double output_value=(*pdf).FullEvaluate
	(base_input, &output_err, gradient);

      //we have weight, error, time to push_back the new event!
      result.push_back(PDF::MC_Event(base_output, output_err, output_value));
      result.back().SetExtra(gradient);

    }

  //now we have to rescale all the weight so they all sum to 1
  //first get the sum of all weight
  double sum_of_weight=0;
  for(unsigned int i=0; i<result.size(); ++i)
    sum_of_weight+=result[i].weight;

  //see if the sum of weights are zero (something wrong??)
  if(sum_of_weight==0)
    cout<<"ERROR: pdf_values from MC_Generator are all zero, results may not make sense"<<endl;
  
  if(sum_of_weight>0)
    {
      //get 1/sum first to make calculation more efficient
      sum_of_weight=1.0/sum_of_weight;
      
      //get gradient scaling factor
      double gradient_rescale=nevt*sum_of_weight/gradient_volume;

      for(unsigned int i=0; i<result.size(); ++i)
	{
	  result[i].weight*=sum_of_weight;
	  result[i].err*=sum_of_weight;
	  result[i].GetExtra()*=gradient_rescale;
	}
    }

  //and that's it! return the list of events
  return result;
}


//generate event quickly
vector<PDF::MC_Event> PDF::MC_Generator::GenerateFast
(const vector<double>& input, int nevt) const
{
  //the number of input must match the input_size
  assert(input.size()==input_size);

  //make sure the PDF MC_Generator refers to exists
  assert(pdf!=NULL);

  //initialize result
  vector<PDF::MC_Event> result;
  result.reserve(nevt);

  //get a list of indices for MC generation
  vector<int> list_of_indices;
  vector<vector<double> > list_of_output;

  list_of_indices.reserve(nevt);
  
  //first get the initial index the input corresponds to
  int initial_index=0;

  int input_index=0;
  for(unsigned int i=0; i<flag.size(); ++i)
    if(flag[i]==INPUT)
      {
	//compute bin_coordinate
	double temp_coordinate=(input[input_index]-pdf->min[i])*
	  pdf->one_over_width[i];

	//compute which bin this falls in
	int temp_bin=pdf->GetBin_helper(input[input_index], i);

	//add the corrdinate index
	initial_index+= temp_bin* (pdf->data.index_helper[i]);
	++input_index;
      }
  list_of_indices.push_back(initial_index);
  
  //now populate all the indices
  for(int i=flag.size()-1; i>=0; --i)
    if(flag[i]==OUTPUT)
      {
	//for all indices already in the list
	//add the possible index + n*helper
	int current_size=list_of_indices.size();
	for(int k=0; k<(pdf->nbins[i]); ++k)
	  for(int j=0; j<current_size; ++j)
	    {
	      list_of_indices.push_back(list_of_indices[j]+ 
					k*pdf->data.index_helper[i]);
	    }
	
      }

  //get a list of nonzero indices
  vector<int> list_of_indices_nonzero;  

  //now grab a bunch of output vectors
  for(unsigned int i=0; i<list_of_indices.size(); ++i)
    {
      //output list does not contain the input
      vector<double> temp_output(pdf->Rank() - input_size);
      vector<int> bin_vec=pdf->data.IndexVec(list_of_indices[i]);
     
      int temp_output_index=0;
      for(unsigned j=0; j<bin_vec.size(); ++j)
	{
	  if(flag[j]==OUTPUT)
	    temp_output[temp_output_index]=
	      pdf->GetBinCenter(bin_vec[j],j);
	  ++temp_output_index;
	}

      //check and see if pdf is nonzero
      
      if(pdf->data[list_of_indices[i]] != 0 || 
	 pdf->data_bias[list_of_indices[i]] != 0)
	{
	  list_of_output.push_back(temp_output);
	  list_of_indices_nonzero.push_back(list_of_indices[i]);
	}
    }

  //now we got a full list of indices
  //time to do some MC!
  
  for(int i=0; i<nevt; ++i)
    {
      //generate an int
      int random_int = irand() % list_of_indices_nonzero.size();
      
      //push back event and that's it!
      result.push_back
	(PDF::MC_Event(list_of_output[random_int], 
		       pdf->data_err[list_of_indices_nonzero[random_int]],
		       pdf->data[list_of_indices_nonzero[random_int]]));
    }


  //now we have to rescale all the weight so they all sum to 1
  //first get the sum of all weight
  double sum_of_weight=0;
  for(unsigned int i=0; i<result.size(); ++i)
    sum_of_weight+=result[i].weight;

  //see if the sum of weights are zero (something wrong??)
  if(sum_of_weight==0)
    cout<<"ERROR: pdf_values from MC_Generator are all zero, results may not make sense"<<endl;
  
  if(sum_of_weight>0)
    {
      //get 1/sum first to make calculation more efficient
      sum_of_weight=1.0/sum_of_weight;
      
      for(unsigned int i=0; i<result.size(); ++i)
	{
	  result[i].weight*=sum_of_weight;
	  result[i].err*=sum_of_weight;
	}
    }
  
  //and that's it
  return result;
}

vector<PDF::MC_Event> PDF::MC_Generator::GenerateFull
(const vector<double>& input) const
{
  //the number of input must match the input_size
  assert(input.size()==input_size);

  //make sure the PDF MC_Generator refers to exists
  assert(pdf!=NULL);

  //initialize result
  vector<PDF::MC_Event> result;

  //get a list of indices for MC generation
  vector<int> list_of_indices;
  vector<vector<double> > list_of_output;

  
  //first get the initial index the input corresponds to
  int initial_index=0;

  int input_index=0;
  for(unsigned int i=0; i<flag.size(); ++i)
    if(flag[i]==INPUT)
      {
	//compute bin_coordinate
	double temp_coordinate=(input[input_index]-pdf->min[i])*
	  pdf->one_over_width[i];

	//compute which bin this falls in
	int temp_bin=pdf->GetBin_helper(input[input_index], i);

	//add the corrdinate index
	initial_index+= temp_bin* (pdf->data.index_helper[i]);
	++input_index;

      }
  list_of_indices.push_back(initial_index);
  
  //now populate all the indices
  for(int i=flag.size()-1; i>=0; --i)
    if(flag[i]==OUTPUT)
      {
	//for all indices already in the list
	//add the possible index + n*helper
	int current_size=list_of_indices.size();
	for(int k=0; k<(pdf->nbins[i]); ++k)
	  for(int j=0; j<current_size; ++j)
	    {
	      list_of_indices.push_back(list_of_indices[j]+ 
					k*pdf->data.index_helper[i]);
	    }

      }

  //now grab a bunch of output vectors
  for(unsigned int i=0; i<list_of_indices.size(); ++i)
    {
      //output list does not contain the input
      vector<double> temp_output(pdf->Rank() - input_size);
      vector<int> bin_vec=pdf->data.IndexVec(list_of_indices[i]);
     
      int temp_output_index=0;
      for(unsigned j=0; j<bin_vec.size(); ++j)
	{
	  if(flag[j]==OUTPUT)
	    temp_output[temp_output_index]=
	      pdf->GetBinCenter(bin_vec[j],j);
	  ++temp_output_index;
	}
      list_of_output.push_back(temp_output);
    }

  //now we got a full list of indices
  //since this is GenerateFull, we simply dump all the values out
  
  for(int i=0; i<list_of_output.size(); ++i)
    {     
      //push back event and that's it!
      result.push_back
	(PDF::MC_Event(list_of_output[i], 
		       pdf->data_err[list_of_indices[i]],
		       pdf->data[list_of_indices[i]]));
    }


  //rescale the weights if desired
  /*
  
  //now we have to rescale all the weight so they all sum to 1
  //first get the sum of all weight
  double sum_of_weight=0;
  for(unsigned int i=0; i<result.size(); ++i)
    sum_of_weight+=result[i].weight;

  //see if the sum of weights are zero (something wrong??)
  if(sum_of_weight==0)
    cout<<"ERROR: pdf_values from MC_Generator are all zero, results may not make sense"<<endl;
  
  if(sum_of_weight>0)
    {
      //get 1/sum first to make calculation more efficient
      sum_of_weight=1.0/sum_of_weight;
      
      for(unsigned int i=0; i<result.size(); ++i)
	{
	  result[i].weight*=sum_of_weight;
	  result[i].err*=sum_of_weight;
	}
    }
  
  */

  //and that's it
  return result;
}



PDF::MC_Generator PDF::Generator(const vector<MC_Flag>& flag_) const
{
  return PDF::MC_Generator(flag_, this);
}
